/******************************************************//**
 * C++ library of the Lévy-Gradient Descent (Lévy-GD) method.
 *
 * Copyright (c) 2020-2031 Yi Zhang (zhangyiss@icloud.com)
 * All rights reserved.
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in
 * all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 * THE SOFTWARE.
 *********************************************************/

#ifndef _LGD_H
#define _LGD_H

/**
 * @brief    A simple definition of the float type we use here. 
 * Easy to change in the future. For now it is just an alias of the double type.
 */
typedef double lgd_float;

/**
 * @brief    Parameters of the L-GD method.
 */
struct lgd_para
{
	/**
	 * Maximal times of the lévy flight. The default value is 1000. The user can 
	 * adjust this parameter by passing a lgd_para type to the lgd_solver() 
	 * functions etc.
	 */
	int flight_times;

	/**
	 * Epsilon for convergence test. This parameter determines the accuracy 
	 * with which the solution is to be found. If this parameter is negative. 
	 * the iteration process will stop till the maximal flight times is reached. 
	 * The default is -1.0.
	 */
	lgd_float epsilon;

	/**
	 * Standard deviation of v that is used to calculate the distance of each
	 * lévy flight in length = stddev_u/|stddev_v|^{1/beta}. This parameter is 
	 * typically given as 1.0.
	 */
	lgd_float stddev_v;

	/**
	 * Scale parameter for calculating stddev_u and the flying length. Must be at 
	 * (1.0, 2.0). The default value is 1.5. The bigger beta is the smaller of the 
	 * range of flying length gets.
	 */
	lgd_float beta;

	/**
	 * scale parameter multiplied by the flying length. The default value is 0.01. 
	 * The parameter should be set according to the expected convergence speed. Normally, 
	 * The bigger alpha is, the faster the L-GD convergences. However, the L-GD may 
	 * miss the optimized solutions if alpha was too big.
	 */
	lgd_float alpha;

	/**
	 * Sigma for the stagnation point test. The algorithm will take one search 
	 * orthogonal with the last iteration if the module of the gradients is smaller 
	 * than sigma. This mechanism helps the algorithm escaping from stagnation 
	 * points such as local minimal or saddle points.The default is 1e-8.
	 */
	lgd_float sigma;

	/**
	 * module length of the search vectors for calculating pseudo gradients. 
	 * The default value is 1.0.The parameter should be set according to the scale 
	 * of the solution.
	 */
	lgd_float mu;

	/**
	 * number of the orthogonal vectors used for calculating pseudo gradients. 
	 * The default value equals to the size of solution.The parameter must be 
	 * equal or smaller than the size of the solution. The minimal value is 2.
	 */
	int orth_num;

	/**
	 * Whether to calculate pseudo model gradients using the in-build subroutine. It takes 
	 * 2*orth_num times forward evaluations to calculate the gradients for every iteration. 
	 * Therefore, it may cost significantly more time to search for the optimized solution. 
	 * If the parameter is set to false, then the model gradients must be provided by the 
	 * evaluation function. The default value is 0. 
	 */
	int pseudo_grad;
};

/**
 * @brief    Callback interface for calculating the value of objective function
 * and the corresponding model gradients.
 * 
 * @note     The users should define their own evaluation functions for specified optimization 
 * problems. Once the problem is defined, this function could be passed to the lgd_solver() or 
 * slgd_solver() functions for evaluations.
 * 
 * @param    instance   The user data sent for the lgd_solver() and slgd_solver() functions. If the 
 * evaluation function is a member function of a class, this variable should be set to the 'this' pointer.
 * @param    x          Pointer of the solution array.
 * @param    g          Pointer of the model gradient array. This array is created by the lgd_solver() or 
 * slgd_solver() function. So don't free the memory yourself.
 * @param    n_size     Length of the solution.
 * @param    m          Index of the observation.
 * 
 * @return   Value of objective function.
 */
typedef lgd_float (*lgd_evaulate_ptr)(void *instance, const lgd_float *x, lgd_float *g, 
	const int n_size, const int m);

/**
 * @brief    Callback interface for monitoring the progress and terminate the iteration 
 * if it is necessary.
 * 
 * @note     The users could define their own monitoring functions for specified optimization 
 * problems. Once the function is defined, it could be passed to the lgd_solver() or 
 * slgd_solver() functions. Passing a NULL pointer to the solver functions will simply disable 
 * the monitoring feature.
 * 
 * @param    instance   The user data sent for the lgd_solver() and slgd_solver() functions. If the 
 * evaluation function is a member function of a class, this variable should be set to the 'this' pointer.
 * @param    fx         Value of objective function currently recorded.
 * @param    best_x     Best solution found currently. This pointer will be NULL for the slgd_solver() function.
 * @param    param      Parameters applied to the L-GD algorithm.
 * @param    n_size     Length of the solution array.
 * @param    curr_k     Times of the current flight.
 * @param    best_k     Search times of which the best_fx is recored.
 * 
 * @return   int        Zero to continue the searching process. Otherwise, the optimization 
 * process will be terminated.
 */
typedef int (*lgd_progress_ptr)(void *instance, lgd_float fx, const lgd_float *best_x, 
	const lgd_para *param, const int n_size, const int curr_k, const int best_k);

/**
 * @brief      Locate memory for a ldg_float pointer type.
 *
 * @param[in]  n_size  Size of the ldg_float array.
 *
 * @return     Pointer of the data
 */
lgd_float *lgd_malloc(const int n_size);

/**
 * @brief      Destroy memory used by the lgd_float type array.
 *
 * @param      x     Pointer of the array.
 */
void lgd_free(lgd_float *x);

/**
 * @brief      Return a lgd_para type instance with default values.
 *
 * @return     A lgd_para type instance.
 */
lgd_para lgd_default_parameters();

/**
 * @brief      Display or throw out a string explanation for the solver functions' return values.
 *
 * @param[in]  er_index  The error index returned by the solver functions.
 * @param[in]  er_throw  Throw out a string of the explanation.
 */
void lgd_error_str(int er_index, bool er_throw = false);

/**
 * @brief      Calculate step lengths of a lévy flight
 * 
 * @warning    You must create the array before you call this function. And the 
 * size of the array must be equal to the flight_times.
 *
 * @param[in]  param  Pointer of a lgd_para object
 * @param      x      The return of a lgd_float array with the size of flight_times
 *
 * @return     Status of the function
 */
int lgd_levy_length(const lgd_para *param, lgd_float *x);

/**
 * @brief      A Lévy-Gradient Descent (LGD) solver function.
 * 
 * @note       The size of all arrays must be equal to n_size. Expect the size of trace is n_size*flight_times.
 *
 * @param      trace         Record solutions generated by the search strategy. Use nullptr to disable the feature. This array
 * should be an n_size array of the lgd_float pointer type. Any pointers that are not nullptr should be initialized a lgd_float 
 * array of the sizes equal to flight_times. And temporary solutions of model parameters at these indexes will be stored.
 * @param      m_mean        Averaged solution calculated using all previous solutions. Use nullptr to disable the feature.
 * @param      m_stddev      Standard deviations of the averaged solutions. Use nullptr to disable the feature.
 * @param[in]  Evafp         Callback function for calculating the objective function and its gradient.
 * @param[in]  Profp         Callback function for monitoring the optimization process.
 * @param      fx            Returned best value of the objective function by now.
 * @param      m             Pointer of the solution array. The best found solution will be returned by this pointer.
 * @param[in]  n_size        Length of the solution array.
 * @param[in]  param         Parameters of L-GD algorithm.
 * @param      instance      The user data sent for the function by the client. If the solver function 
 * is a member function of a class, this variable should be set to the 'this' pointer.
 * @param      low_limit     Lower boundary of the solution. The default is nullptr.
 * @param      high_limit    Upper boundary of the solution. The default is nullptr.
 * @param      alpha_range   Absolute values of the estimated changing range of the solution. The default assumes
 * that the ranges of the variables of the solution are all around 1.0. This array will be calculated automatically 
 * if the upper and lower bounds of the solution are given. Given the estimated ranges will overwrite the array.
 *
 * @return     Status of the function.
 */
int lgd_solver(lgd_float **trace, lgd_float *m_mean, lgd_float *m_stddev, lgd_evaulate_ptr Evafp, 
	lgd_progress_ptr Profp, lgd_float *fx, lgd_float *m, const int n_size, const lgd_para *param, 
	void *instance, const lgd_float *low_limit = nullptr, const lgd_float *high_limit = nullptr, 
	const lgd_float *alpha_range = nullptr);

/**
 * @brief      A Stochastic Lévy-Gradient Descent (SLGD) solver function.
 * 
 * @note       The size of all arrays must be equal to n_size. Expect the size of trace is n_size*m_size*flight_times.
 *
 * @param      trace         Record solutions generated by the search strategy. Use nullptr to disable the feature. This array
 * should be an n_size array of the lgd_float pointer type. Any pointers that are not nullptr should be initialized a lgd_float 
 * array of the sizes equal to flight_times*m_size. And temporary solutions of model parameters at these indexes will be stored.
 * @param      m_mean        Averaged solution calculated using solutions from the last flight.
 * @param      m_stddev      Standard deviations of the averaged solutions.
 * @param[in]  Evafp         Callback function for calculating the objective function and its gradient.
 * @param[in]  Profp         Callback function for monitoring the optimization process.
 * @param      fx            Returned best value of the objective function by now.
 * @param      m             Pointer of the solution array. The best found solution will be returned by this pointer.
 * @param[in]  n_size        Length of the solution array.
 * @param[in]  m_size        Length of the observation array.
 * @param[in]  param         Parameters of L-GD algorithm.
 * @param      instance      The user data sent for the function by the client. If the solver function 
 * is a member function of a class, this variable should be set to the 'this' pointer.
 * @param      low_limit     Lower boundary of the solution. The default is nullptr.
 * @param      high_limit    Upper boundary of the solution. The default is nullptr.
 * @param      alpha_range   Absolute values of the estimated changing range of the solution. The default assumes
 * that the ranges of the variables of the solution are all around 1.0. This array will be calculated automatically 
 * if the upper and lower bounds of the solution are given. Given the estimated ranges will overwrite the array.
 *
 * @return     Status of the function.
 */
int slgd_solver(lgd_float **trace, lgd_float *m_mean, lgd_float *m_stddev, lgd_evaulate_ptr Evafp, 
	lgd_progress_ptr Profp, lgd_float *fx, lgd_float *m, const int n_size, const int m_size, 
	const lgd_para *param, void *instance, const lgd_float *low_limit = nullptr, 
	const lgd_float *high_limit = nullptr, const lgd_float *alpha_range = nullptr);

/**
 * @brief        A virtual class interface for the L-GD algorithm
 * 
 * @warning      You should not edit this class manually.
 * 
 * To use this virtual solver class, you must define your own solver class for any specified 
 * optimization problems and inherit this virtual class as a father class. Moreover, you need to 
 * defined the Evaluate() function for the optimization question. Subsequently, the Minimize() or 
 * MinimizeStochastic() function could be called to search for the optimized solutions. Parameters 
 * of the L-GD algorithm could be adjusted by passing a lgd_para object to the set_lgd_parameter() 
 * function. And customized monitoring functions could be defined by overloading the 
 * Progress() function.
 */
class LGD_Solver
{
protected:
	lgd_para param_;

public:
	LGD_Solver();
	virtual ~LGD_Solver(){}

	static lgd_float _Evaluate(void *instance, const lgd_float *x, lgd_float *g, 
		const int n_size, const int m)
	{
		return reinterpret_cast<LGD_Solver*>(instance)->Evaluate(x, g, n_size, m);
	}

	/**
	 * @brief      Interface of the evaluation function.
	 * 
	 * @note       This virtual function only declares the interface of the evaluation function. The user 
	 * must define this function to use the minimize functions.
	 *
	 * @param[in]  x       Pointer of the solution array.
	 * @param      g       Pointer of the model gradient array. This array is created by the lgd_solver() or 
	 * slgd_solver() function. So don't free the memory yourself.
	 * @param[in]  n_size  Length of the solution.
	 * @param[in]  m       Index of the observation.
	 *
	 * @return     Value of objective function.
	 */
	virtual lgd_float Evaluate(const lgd_float *x, lgd_float *g, const int n_size, const int m) = 0;

	static int _Progress(void *instance, lgd_float fx, const lgd_float *best_x, 
		const lgd_para *param, const int n_size, const int curr_k, const int best_k)
	{
		return reinterpret_cast<LGD_Solver*>(instance)->Progress(fx, best_x, param, n_size, curr_k, best_k);
	}

	/**
	 * @brief      Monitoring function of the L-GD algorithm.
	 *
	 * @param[in]  fx      Value of objective function currently recorded.
	 * @param[in]  best_x  Best solution found currently. This pointer will be NULL for the slgd_solver() function.
	 * @param[in]  param   Parameters applied to the L-GD algorithm.
	 * @param[in]  n_size  Length of the solution array.
	 * @param[in]  curr_k  Times of the current flight.
	 * @param[in]  best_k  Search times of which the best_fx is recored.
	 *
	 * @return     Zero to continue the searching process. Otherwise, the optimization 
	 * process will be terminated.
	 */
	virtual int Progress(lgd_float fx, const lgd_float *best_x, 
		const lgd_para *param, const int n_size, const int curr_k, const int best_k);

	/**
	 * @brief      Sets parameters of the LGD algorithm using an external variable
	 *
	 * @param[in]  in_param  In put parameter
	 */
	void set_lgd_parameter(const lgd_para &in_param);

	/**
	 * @brief      The L-GD solver function.
	 *
	 * @param      trace         Record solutions generated by the search strategy. Use nullptr to disable the feature. This array
	 * should be an n_size array of the lgd_float pointer type. Any pointers that are not nullptr should be initialized a lgd_float 
	 * array of the sizes equal to flight_times. And temporary solutions of model parameters at these indexes will be stored.
	 * @param      m_mean        Averaged solution calculated using all previous solutions. Use nullptr to disable the feature.
	 * @param      m_stddev      Standard deviations of the averaged solutions. Use nullptr to disable the feature.
	 * @param      fx            Returned best value of the objective function by now.
	 * @param      m             Pointer of the solution array. The best found solution will be returned by this pointer.
	 * @param[in]  n_size        Length of the solution array.
	 * @param[in]  low_limit     Lower boundary of the solution. The default is nullptr.
	 * @param      high_limit    Upper boundary of the solution. The default is nullptr.
	 * @param      alpha_range   Absolute values of the estimated changing range of the solution. The default assumes
	 * that the ranges of the variables of the solution are all around 1.0. This array will be calculated automatically 
	 * if the upper and lower bounds of the solution are given. Given the estimated ranges will overwrite the array.
	 * @param[in]  verbose       Whether to set the function to verbose mode to receive all return messages. If this variable is 
	 * set to false. The function will only display messages when something goes wrong.
	 * @param[in]  er_throw      Whether to throw out a error message instead of displaying it on the screen.
	 */
	void Minimize(lgd_float **trace, lgd_float *m_mean, lgd_float *m_stddev, lgd_float *fx, 
		lgd_float *m, const int n_size, const lgd_float *low_limit = nullptr, 
		const lgd_float *high_limit = nullptr, const lgd_float *alpha_range = nullptr, 
		bool verbose = true, bool er_throw = false);

	/**
	 * @brief      The stochastic L-GD solver function.
	 *
	 * @param      trace         Record solutions generated by the search strategy. Use nullptr to disable the feature. This array
	 * should be an n_size array of the lgd_float pointer type. Any pointers that are not nullptr should be initialized a lgd_float 
	 * array of the sizes equal to flight_times*m_size. And temporary solutions of model parameters at these indexes will be stored.
	 * @param      m_mean        Averaged solution calculated using all previous solutions. Use nullptr to disable the feature.
	 * @param      m_stddev      Standard deviations of the averaged solutions. Use nullptr to disable the feature.
	 * @param      fx            Returned best value of the objective function by now.
	 * @param      m             Pointer of the solution array. The best found solution will be returned by this pointer.
	 * @param[in]  n_size        Length of the solution array.
	 * @param[in]  m_size        Length of the observation array.
	 * @param[in]  low_limit     Lower boundary of the solution. The default is nullptr.
	 * @param      high_limit    Upper boundary of the solution. The default is nullptr.
	 * @param      alpha_range   Absolute values of the estimated changing range of the solution. The default assumes
	 * that the ranges of the variables of the solution are all around 1.0. This array will be calculated automatically 
	 * if the upper and lower bounds of the solution are given. Given the estimated ranges will overwrite the array.
	 * @param[in]  verbose       Whether to set the function to verbose mode to receive all return messages. If this variable is 
	 * set to false. The function will only display messages when something goes wrong.
	 * @param[in]  er_throw      Whether to throw out a error message instead of displaying it on the screen.
	 */
	void MinimizeStochastic(lgd_float **trace, lgd_float *m_mean, lgd_float *m_stddev, 
		lgd_float *fx, lgd_float *m, const int n_size, const int m_size, const lgd_float *low_limit = nullptr, 
		const lgd_float *high_limit = nullptr, const lgd_float *alpha_range = nullptr, 
		bool verbose = true, bool er_throw = false);
};

#endif // _LGD_H